%% Relative Coverage Image Analysis Script
% Cassady Rupert, Brown University, 2019-07-05
% Asks user to select a brightfield image for analysis and produces a
% tiled image of binarized subimages and a relative coverage value

%% Select and read in image
% Prompt user to choose image
% Default folder for images is current matlab directory
[file,pathname] = uigetfile('*');
% Assign image filename including full directory path
filename = strcat(pathname,'\',file);
% Read in image file and convert to grayscale
I = imread(filename);
Gscale = rgb2gray(I);

%% Ask user to pick subimage rectangle
% show grayscale image
figure
imshow(Gscale)
title('select ROI that excludes areas outside the well')
% prompt user to select rectangle to analyze
well_rect = getrect;
%% Partition ROI into 16 rectangles to analyze individually. 
% This accounts for changes in brightness due to edge effects from the well
% border and changes in lighting across the well

% crop identified area to have width and height divisible by 4, so
% subimages can be evenly divided to a 4x4 array
width_rem = rem(well_rect(3),4);
if width_rem <1
    width = well_rect(3)-4;
else
    width = well_rect(3)-width_rem-4;
end
height_rem = rem(well_rect(4),4);
if height_rem <1
    height = well_rect(4)-4;
else
    height = well_rect(4)-height_rem-4;
end

%Round identifying corner of drawn rectangle so that it's within the
%selected rectangle
for i=1:2
    well_rect(i) = round(well_rect(i));
end

% Define first subimage of the selected rectangle 
sub_rect = [well_rect(1) well_rect(2) (width/4) (height/4)];

%% Threshold grayscale subimages 
% loop throgh subimages, thresholing each by global method
% Initialize vectors to collect binarized image areas and relative coverage
Area_collect = zeros(16,1);
RelativeCovg_collect = zeros(16,1);
% Initialize counter for collecting subimage values
n=1;
% Initialize figure to collect binzarized subimages
figure
comb_thresh = tight_subplot(4,4,[0 .01],[.01 .01],[.01 .01]);

% Loop through subimages to binarize 16 equally sized sub-rectangles
for j = 1:4
    for i = 1:4
        % Define subimage area in grayscale image
        subI = Gscale(sub_rect(2)+((i-1)*sub_rect(4)):(sub_rect(2)+(i*sub_rect(4))),... 
            sub_rect(1)+((j-1)*sub_rect(3)):(sub_rect(1)+(j*sub_rect(3))));
        % Threshold subimage
        [level,EM] = graythresh(subI);
        BW = imbinarize(subI,level);
        % Calculate and collect subimage area and fraction coverage
        Area_collect(n) = bwarea(BW);
        Sz = size(subI);
        RelativeCovg_collect(n)=Area_collect(n)/(Sz(1)*Sz(2));
        % subplot fills up flipped, so reassign n to m values
        m=j+((i-1)*4);
        axes(comb_thresh(m));imshow(BW);
        n=n+1;
    end
end

% Calculate average image Realtive Coverage
RelativeCovg = mean(RelativeCovg_collect);

% Add annotation to binarized subplot with relative coverage
str = strcat('Relative coverage: ', num2str(RelativeCovg));
annotation('textbox','String',str,'FitBoxToText',...
    'on','BackgroundColor','w');

% Function to customize subplot of binarized subimages
% Adapted from    
% Pekka Kumpulainen 21.5.2012   @tut.fi
% Tampere University of Technology / Automation Science and Engineering
function [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)

% tight_subplot creates "subplot" axes with adjustable gaps and margins
%
% [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
%
%   in:  Nh      number of axes in hight (vertical direction)
%        Nw      number of axes in width (horizontaldirection)
%        gap     gaps between the axes in normalized units (0...1)
%                   or [gap_h gap_w] for different gaps in height and width 
%        marg_h  margins in height in normalized units (0...1)
%                   or [lower upper] for different lower and upper margins 
%        marg_w  margins in width in normalized units (0...1)
%                   or [left right] for different left and right margins 
%
%  out:  ha     array of handles of the axes objects
%                   starting from upper left corner, going row-wise as in
%                   subplot
%        pos    positions of the axes objects
%
%  Example: ha = tight_subplot(3,2,[.01 .03],[.1 .01],[.01 .01])
%           for ii = 1:6; axes(ha(ii)); plot(randn(10,ii)); end
%           set(ha(1:4),'XTickLabel',''); set(ha,'YTickLabel','')

% Pekka Kumpulainen 21.5.2012   @tut.fi
% Tampere University of Technology / Automation Science and Engineering


if nargin<3; gap = .02; end
if nargin<4 || isempty(marg_h); marg_h = .05; end
if nargin<5; marg_w = .05; end

if numel(gap)==1; 
    gap = [gap gap];
end
if numel(marg_w)==1; 
    marg_w = [marg_w marg_w];
end
if numel(marg_h)==1; 
    marg_h = [marg_h marg_h];
end

axh = (1-sum(marg_h)-(Nh-1)*gap(1))/Nh; 
axw = (1-sum(marg_w)-(Nw-1)*gap(2))/Nw;

py = 1-marg_h(2)-axh; 

% ha = zeros(Nh*Nw,1);
ii = 0;
for ih = 1:Nh
    px = marg_w(1);
    
    for ix = 1:Nw
        ii = ii+1;
        ha(ii) = axes('Units','normalized', ...
            'Position',[px py axw axh], ...
            'XTickLabel','', ...
            'YTickLabel','');
        px = px+axw+gap(2);
    end
    py = py-axh-gap(1);
end
if nargout > 1
    pos = get(ha,'Position');
end
ha = ha(:);
end

